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A  pore  scale  model  of  a  polymer  electrolyte  membrane  (PEM)  fuel  cell  cathode  catalyst  layer  is  developed 
which  accounts  for  species  transport,  electrochemical  reactions  and  thermal  transport.  Effective  transport 
parameters  are  computed  over  a  range  of  operating  conditions  including  the  effective  oxygen  diffusivity, 
effective  water  vapor  diffusivity,  effective  proton  conductivity,  effective  electron  conductivity  and  the 
effective  thermal  conductivity.  In  addition,  the  total  amount  of  oxygen  consumption  is  computed  for 
different  operating  conditions.  Finally,  a  critical  assessment  of  the  impact  of  assumptions  made  in  the 
absence  of  detailed  morphological  data  is  presented. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Catalyst  layers  are  at  the  heart  of  PEM  fuel  cells.  They  pro¬ 
vide  pathways  for  the  transport  of  reactants,  protons,  electrons 
and  producs  while  facilitating  oxygen  reduction  at  the  cathode  and 
hydrogen  oxidation  at  the  anode.  Catalyst  layers  are  composed  of 
four  phases:  carbon,  which  allows  for  conduction  of  electrons  and 
support  the  platinum  nanoparticles;  ionomer  binding,  which  pro¬ 
vides  a  path  for  proton  conduction;  platinum  nanoparticles,  which 
provide  sites  for  electrochemical  reactions;  and  pores,  which  allow 
reactant  and  product  gases  to  diffuse  through  the  catalyst  layer. 
All  of  these  phases  are  present  in  an  extremely  thin  layer  (of  the 
order  of  10  pum),  which  is  bounded  on  one  side  by  the  gas  diffusion 
layer  or  microporous  layer  and  on  the  other  side  by  the  polymer 
electrolyte  membrane. 

The  primary  issues  which  prevent  the  commercialization  of  PEM 
fuel  cells  are  high  cost,  degradation  and  poor  cell  performance.  The 
high  cost  is  related  to  the  amount  of  platinum  that  is  used  in  cat¬ 
alyst  layers.  Improving  the  utilization  of  the  platinum  particles  in 
catalyst  layers,  particularly  on  the  cathode  side,  has  the  potential 
to  substantially  reduce  the  cost  of  PEM  fuel  cells  [1].  Degradation 
of  the  ionomer,  carbon  supports  and  platinum  particles  is  also  a 
significant  issue  [2,3].  Understanding  the  mechanisms  by  which 
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these  degradations  occur  can  lead  to  mitigation  strategies  that  can 
improve  the  operating  lifetime  of  fuel  cells.  Finally,  the  poor  cell 
performance  that  is  experienced  by  PEM  fuel  cells  is  due  to  a  large 
activation  polarization  and  mass-transport  losses  at  the  cathode. 
The  activation  polarization  is  caused  by  sluggish  kinetics  of  the 
oxygen  reduction  reaction,  and  is  particularly  large  because  of  the 
low  operating  temperature  of  PEM  fuel  cells.  The  most  severe  mass 
transport  losses  are  caused  by  the  accumulation  of  liquid  water  at 
reaction  sites,  which  is  often  referred  to  as  “flooding”.  Additionally, 
significant  ohmic  losses  through  the  membrane  can  occur  at  dry 
conditions  [4]. 

Because  the  catalyst  layer  is  so  thin,  it  is  very  difficult  to  do 
experimental  studies  which  could  provide  information  that  would 
allow  one  to  reduce  the  cost  and  degradation  or  improve  the  per¬ 
formance  of  a  PEM  fuel  cell.  Because  of  this,  there  has  been  a  surge 
of  interest  in  computational  simulation  and  reconstruction  of  the 
catalyst  layer  at  the  mesoscale  level.  It  is  hoped  that  simulations 
will  eventually  provide  a  cost-effective  approach  for  designing 
optimum  catalyst  layer  microstructures.  However,  there  are  still 
substantial  gaps  in  the  existing  body  of  knowledge  that  prevent 
current  simulation  tools  from  having  these  predictive  capabilities. 
Nonetheless,  mesoscale  simulations  can  provide  valuable  informa¬ 
tion  such  as  the  effective  transport  properties  through  the  catalyst 
layer.  These  properties  could  be  provided  to  a  macroscale  model 
as  part  of  a  multiscale  fuel  cell  model.  In  addition,  important  ideas 
about  multiphase  transport  through  nanoscale  porous  media  could 
be  obtained  from  these  models. 
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The  first  published  work  on  catalyst  layer  simulations  was  done 
by  Wang  [5]  in  two  dimensions.  This  model,  which  included  pro¬ 
ton  transport,  oxygen  transport  and  electrochemical  reactions, 
was  later  extended  to  three  dimensions  and  applied  to  a  regular 
microstructure  [6]  and  a  random  microstructure  [7].  In  the  three- 
dimensional  models,  electron  transport  was  neglected  because 
ion  transport  was  assumed  to  be  the  limiting  factor  for  the  oxy¬ 
gen  reduction  reaction.  The  production  and  transport  of  water 
was  also  accounted  for  in  the  three-dimensional  model.  Addi¬ 
tional  simulations  were  done  on  a  catalyst  layer  structure  that  was 
reconstructed  to  match  a  physical  catalyst  layer  using  stochastic 
methods  [8].  This  work  was  further  extended  to  identify  an  opti¬ 
mum  design  [9],  based  on  a  series  of  parametric  studies  rather  than 
multi-parameter  numerical  optimization  [10,1  ].  These  simulations 
modeled  the  entire  length  of  the  catalyst  layer  and  had  a  resolution 
of  250  nm. 

Kim  and  Pitsch  [11]  used  a  different  approach  to  catalyst  layer 
reconstruction  and  simulation.  They  represented  the  catalyst  layer 
as  consisting  of  carbon  spheres  surrounded  by  an  ionomer  film.  The 
catalyst  layer  reconstruction  was  done  using  a  simulated  annealing 
technique,  where  the  spheres  were  randomly  moved  in  order  to 
obtain  a  prescribed  porosity  and  two-point  autocorrelation  func¬ 
tion.  The  catalyst  layer  was  modeled  with  a  resolution  of  8nm 
and  a  1  pun  x  1  pun  x  1  p.m  section  was  considered.  They  used  a 
Lattice-Boltzmann  method  to  simulate  the  transport  of  oxygen 
and  protons  through  the  catalyst  layer,  but  their  simulation  did 
not  account  for  electrochemical  reactions.  Knudsen  diffusion  was 
accounted  for  by  computing  an  average  pore  radius  which  was 
used  to  compute  a  Knudsen  diffusion  coefficient  for  all  pores  in  the 
solution  domain.  They  obtained  effective  transport  parameters  for 
oxygen  and  protons,  and  expressed  the  relationship  between  vol¬ 
ume  fractions  and  effective  transport  parameters  by  using  a  power 
law  approximation. 

More  recently,  a  small  section  of  a  PEM  fuel  cell  catalyst  layer 
was  reconstructed  and  modeled  at  a  resolution  of  2nm  [12].  A 
random  reconstruction  algorithm  which  attempted  to  account  for 
agglomeration  was  considered  in  the  model.  In  addition,  a  different 
Knudsen  number  was  used  for  each  pore  cell  based  on  the  diame¬ 
ters  in  the  x-,  y-  and  z-directions.  Effective  transport  properties  for 
oxygen  diffusivity  and  proton  conductivity  were  computed  for  five 
different  geometries.  Specific  attention  was  focused  on  the  number 
of  seeds  available  for  agglomerates  to  grow  on. 

In  a  previous  paper  [13],  the  authors  used  a  stochastic  approach 
to  reconstruct  a  section  of  a  PEM  fuel  cell  catalyst  layer  and  mod¬ 
eled  oxygen  diffusion,  proton  conduction,  electron  conduction  and 
electrochemical  reactions  for  a  wide  range  of  random  microstruc¬ 
tures.  Hundreds  of  different  microstructures  were  used,  resulting  in 
a  large  number  of  effective  transport  parameters  at  different  pore, 
ionomer  and  carbon  volume  fractions.  The  microstructural  parame¬ 
ters  for  the  reconstructed  catalyst  layer  were  changed  to  determine 
how  this  would  affect  effective  transport  parameters  and  total  con¬ 
sumption  values.  One  of  the  key  findings  of  this  work  was  that  when 
larger  carbon  spheres  are  used  in  the  catalyst  layer,  higher  effective 
oxygen  diffusivities  and  proton  conductivities  are  obtained. 

The  current  work  builds  on  the  previous  work  [13]  and  extends 
the  model  to  account  for  water  vapor  concentration  and  temper¬ 
ature.  The  effects  of  water  vapor  diffusion,  electro-osmotic  drag, 
thermal  conduction,  ohmic  heating,  and  activation  losses  are  con¬ 
sidered.  The  fully  coupled  set  of  discretized  governing  equations  is 
solved  for  a  wide  range  of  microstructures  at  various  temperatures, 
oxygen  concentrations,  and  water  vapor  concentrations.  Effective 
oxygen  diffusivities,  effective  water  vapor  diffusivities,  effective 
proton  conductivities  and  effective  electron  conductivities  are 
computed  for  a  wide  range  of  microstructures.  In  addition,  effective 
thermal  conductivities  for  different  microstructures  are  computed. 
To  the  knowledge  of  the  authors,  this  is  the  first  mesoscale  PEM 


Fig.  1.  Typical  stochastically  reconstructed  catalyst  layer  section 
(200  nm  x  200  nm  x  200  nm).  Carbon  spheres  are  in  blue,  ionomer  is  green, 
and  pores  are  orange  and  red. 


fuel  cell  catalyst  layer  simulation  to  obtain  effective  thermal  con¬ 
ductivity  data. 

2.  Numerical  method 

2.1.  Catalyst  layer  reconstruction 

A  stochastic  method  is  used  to  reconstruct  random  microstruc¬ 
tures  which  represent  pieces  of  the  catalyst  layer.  Carbon  spheres 
are  randomly  placed  in  the  domain  subject  to  several  constraints. 
Each  carbon  sphere  has  a  probability  that  it  will  be  required  to  over¬ 
lap  with  a  previously  placed  sphere,  but  it  can  only  overlap  by  an 
amount  that  is  less  than  the  specified  overlap  tolerance.  Once  all 
the  carbon  spheres  have  been  placed  in  the  domain,  an  ionomer 
layer  of  a  given  thickness  is  extruded  from  the  sphere  surfaces. 
Thus  the  domain  consists  of  three  phases:  carbon,  ionomer,  and 
pores.  More  details  of  the  reconstruction  algorithm  have  been  given 
in  previous  works  [14,13]  and  are  not  repeated  here.  In  the  cur¬ 
rent  work,  an  ionomer  thickness  of  4  nm  is  used  along  with  carbon 
sphere  radii  of  12nm  and  20  nm.  An  overlap  probability  of  0.9  is 
used  (representing  the  tendency  of  carbon  particles  to  agglomer¬ 
ate)  and  an  overlap  tolerance  of  0.01  is  used  (rigid  carbon  spheres 
will  not  have  a  significant  overlap).  The  domain  size  considered 
in  this  work  is  200  nm  x  200  nm  x  200  nm,  which  is  discretized 
with  a  100x100x100  structured  mesh.  The  stochastic  simula¬ 
tions  performed  in  this  study  required  the  generation  of  hundreds 
of  different  computational  samples;  an  example  is  shown  in  Fig.  1. 

2.2.  Governing  equations 

The  model  incorporates  governing  equations  for  oxygen,  water 
vapor,  protons,  electrons,  and  temperature.  Steady-state  operating 
conditions  are  assumed  and  gas  transport  is  assumed  to  take  place 
through  binary  diffusion.  The  effects  of  liquid  water  are  not  con¬ 
sidered  in  the  model.  The  electrochemical  reactions  are  assumed 
to  take  place  uniformly  at  the  carbon/ionomer  interface  and  the 
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Tafel  approximation  is  used,  which  is  valid  at  the  cathode  for  high 
current  densities. 

The  continuity  equation  for  oxygen  is  composed  of  two 
expressions  which  represent  diffusion  and  consumption  through 
electrochemical  reactions,  and  is  given  as 


V- 


^-Do2VCo2  +  1  int 


1  .  Cq2 

4F,0CO2,re/ 


exp 


f-cccF 
V  RT 


=  0,  (1) 


where  c02  is  the  oxygen  concentration,  F  is  Faraday’s  constant, 
co2,ref  is  the  oxygen  concentration  at  standard  temperature  and 
pressure,  R  is  the  universal  gas  constant  and  T  is  the  temperature. 
The  first  part  of  Eq.  (1 )  represents  oxygen  diffusion.  Because  of  the 
small  pore  diameters  in  the  catalyst  layer,  Knudsen  diffusion  must 
be  accounted  for  in  the  pore  regions.  The  Knudsen  diffusivity  for 
oxygen  is  computed  according  to  the  dusty  gas  model  [15,16]  as 
[17-19] 


£>o2,/<n(cm2s  1 )  =  4850dp 


(2) 


where  dp  is  the  local  pore  diameter  for  a  given  pore  cell,  which  is 
computed  by  taking  an  average  of  13  different  lengths  as  described 
in  a  previous  work  [13].  The  bulk  diffusivity  is  given  according  to 
binary  diffusion  of  oxygen  into  nitrogen  [20]  as 


Do2,b  (cm2  s  1 


=  2.2  x  10-1 


(7/(293.2 1<))3/2 
p/(l  atm) 


(3) 


where  p  is  the  pressure  of  the  catalyst  layer,  which  is  taken  to  be 
2  atm  in  the  current  work.  The  total  oxygen  diffusivity  in  the  pore 
region  is  computed  using  the  Bosanquet  approximation  [21  ]  as 


11  1 

Dq2  Do2,i<n  Do,* 


(4) 


The  carbon  phase  is  assumed  to  be  impermeable  to  oxygen,  but 
the  oxygen  diffusivity  through  Nation  which  has  been  recast  at  a 
temperature  of  120  °C  is  taken  from  experimental  data  [22]  to  be 

Dq2  ,  (cm2  s”1 )  =  10“6  X  (0.1543(7  -  273)  -  1.65).  (5) 


Table  1 

Curve-fit  coefficients  and  parameter  values  used  in  the  simulation. 


Coefficient 

Value 

Coefficient 

Value 

Cl 

2.8133  x  10“4 

Pi 

-2846.4 

c2 

1.328355 

P2 

-411.24 

C3 

-1.1642xl0-2 

P3 

-10.554 

c4 

3.442175  x  10“5 

P4 

0.16636 

C5 

-3.33815  xlO-8 

C6 

-7.2939  x  10“4 

where  cH2o  is  the  water  vapor  concentration,  nd  is  the  electro- 
osmotic  drag  coefficient,  and  crm  is  the  membrane  conductivity.  In 
the  pore  region,  the  expressions  for  the  Knudsen  diffusivity  of  water 
vapor  [17-19]  and  the  bulk  diffusivity  of  water  vapor  [20]  are  given 
as: 


D  H2o,/cn(cm2s  1 )  =  4850dj 


T_ 

18’ 


(10) 


,3/2 


F^H2o,b  (cm2  s  l)  =  2.93x10 


_!  (T/298.2  K) 
p/(l  atm) 


(11) 


The  Bosanquet  approximation  [21  ]  is  used  to  compute  the  total 
water  vapor  diffusivity  in  the  pore  region  as 


DH20  Dh20, I<n  DH20,b 

For  the  ionomer  region,  the  water  vapor  diffusivity  is  given  as 
[24] 

HH2o,i(cm2  s-1)  =  6.0  x  10-6.  (13) 

The  water  vapor  diffusivity  in  the  carbon  region  is  set  to  zero. 
The  value  for  the  electro-osmotic  drag  coefficient  nd  is  taken  to  be 
1  [25]. 

The  proton  continuity  equation  consists  of  two  terms  which 
represent  proton  conduction  and  consumption  of  protons  through 
electrochemical  reactions.  It  is  given  as 


The  second  part  of  Eq.  (1)  refers  to  the  consumption  of  oxy¬ 
gen  through  electrochemical  reactions.  The  term  \  int  is  included  to 
denote  that  this  term  is  only  nonzero  at  the  interfaces  of  carbon 
and  ionomer  cells.  In  this  respect,  the  model  considers  an  ideal¬ 
ized  scenario  where  a  very  thin  platinum  layer  uniformly  covers 
the  carbon  spheres  and  the  electrochemical  reaction  takes  place 
at  every  carbon/ionomer  interface.  As  the  model  is  further  devel¬ 
oped,  discrete  platinum  particles  will  be  considered  along  with  a 
correlation  between  the  platinum  loading  and  the  number  of  plat¬ 
inum  particles.  The  exchange  current  density  i'o  and  the  cathodic 
charge  transfer  coefficient  ac  are  given  according  to  experimental 
data  from  a  platinum/Nafion  electrode  [23]  as 


i0(Acm-2)  =  (0.029825)10“1521'93/t, 

(6) 

ac  =  (2.25  x  10“3)7-  1.78  x  10-1. 

(7) 

The  overpotential  is  given  as 

P  =  <l>s  ~  0m  > 

(8) 

where  0S  is  the  solid  potential  and  </>m  is  the  membrane  potential. 

The  governing  equation  for  water  vapor  includes  terms  for 
diffusion,  electro-osmotic  drag,  and  production  through  electro¬ 
chemical  reactions.  It  is  given  as 

V  ■  (-Dh2oVch2o  -  -  lint 


V- 


CTjtj  V0m  +  1  jnt 


(14) 


The  membrane  conductivity  is  dependent  upon  the  temperature 
and  water  vapor  concentration.  For  this  simulation,  it  is  assumed 
that  the  membrane  is  recast  at  1 20  °C,  and  a  curve-fit  of  experimen¬ 
tal  data  [26]  is  used  to  compute  the  membrane  conductivity  in  the 
ionomer  phase  as 


<7m  (Scm  1 )  =  Ci  exp([c2T  -  c3T2  +  c4T3  -  c5T4]a)  +  c6.  (15) 


The  relative  humidity  a  is  computed  as 
Ch2o  RT 

Psat 


(16) 


while  the  saturation  pressure  Psat  is  computed  as 

Psat  (Pa)  =  PI  +p2(7  -  273)  -  p3(7  -  273 )2  +p4(7  -  273 )3.  (17) 


The  curve-fitting  parameters  for  Eqs.  (15)  and  (17)  are  given  in 
Table  1. 

The  electron  continuity  equation  consists  of  a  conduction  and  a 
consumption  term  and  is  given  as 


V- 


V  <ps  +  1  int 


(18) 


JLi  C°2 

2 F  °C02,ref 


=  0, 


(9) 


The  solid  conductivity  crs  is  10 Scm-1  [27]  in  the  carbon  region 
and  is  zero  in  the  pore  and  ionomer  regions. 
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Table  2 

Thermal  conductivities  for  different  regions  of  the  catalyst  layer. 


Region 

X  (Wcm-1  K1 

) 

Pore 

2.99 

Ionomer 

20.0  [28] 

Carbon 

37.5  [29] 

Table  3 

Test  matrix  of  boundary  conditions. 

Parameter  Set  1 

Set  2 

Boundary  1 

Boundary  2 

Boundary  1 

Boundary  2 

T  (K)  343.0 

343.001 

353.0 

353.001 

Co2  (mol  irr3)  14.1 

14.0 

7.1 

7.0 

Ch2o  (mol  itt3)  10.0 

10.1 

8.0 

8.1 

The  heat  equation  consists  of  a  conduction  term,  ohmic  heating 
terms,  and  reaction  heating  terms  and  is  given  as 


V  •  (— AVT)  = 


(V<fe)2  |  (V(/>m)2 
crs  crm 


-1  f„tV- 


io 


co2 

c0  2,  ref 


exp 


f-dcF 
V  RT 


(rj  +  n), 


(19) 


where  X  is  the  thermal  conductivity.  The  thermal  conductivities 
used  in  each  region  are  shown  in  Table  2. 

The  Peltier  coefficient  77  for  the  oxygen  reduction  reaction  can 
be  approximated  as  [30] 

(2°) 

where  A is  the  entropy  change  of  the  half-cell  reaction,  which  is 
taken  to  be  326.6 J  mol-1  K_1  [31  ]  for  the  simulation. 


2.3.  Boundary  conditions 

In  the  three-dimensional  domain,  Dirichlet  boundary  conditions 
are  applied  at  two  opposite  faces  of  the  cubic  domain,  while  peri¬ 
odic  boundary  conditions  are  imposed  at  the  other  four  faces.  A  test 
matrix  of  simulation  cases  is  considered  with  a  number  of  varying 
boundary  conditions,  which  are  shown  in  Table  3.  A  small  tempera¬ 
ture  difference  of  0.001  K  between  opposite  boundaries  is  imposed, 
while  oxygen  and  water  vapor  concentrations  differ  by  0.1  mol  m-3 
on  opposite  boundaries.  In  addition  to  the  boundary  conditions 
listed  in  Table  3,  cases  are  run  on  microstructures  with  sphere  radii 
of  12  nm  and  microstructures  with  sphere  radii  of  20  nm.  Thus  16 
different  sets  of  simulation  data  are  run  using  different  sphere  radii 
and  the  different  combinations  of  data  considered  in  Table  3. 

The  ion  and  electron  potential  boundary  conditions  are  the  same 
for  every  simulation  case  and  are  given  in  Table  4.  In  the  cathode 
of  a  PEM  fuel  cell  catalyst  layer,  the  proton  and  electron  current  is 
driven  by  the  consumption  of  protons  and  electrons.  The  electrons 
have  no  place  to  go  once  they  reach  the  membrane  side  of  the  cat¬ 
alyst  layer  and  the  protons  cannot  go  anywhere  once  they  reach 
the  GDL  side  of  the  catalyst  layer.  Thus,  their  potentials  were  speci¬ 
fied  such  that  the  proton  and  electron  currents  through  the  random 
microstructures  are  approximately  equal  and  the  overpotential  at 
each  boundary  is  close  to  0.4  V  (typical  for  high  current  densities). 


Table  4 

Potential  boundary  conditions. 


Variable 

Boundary  1 

Boundary  2 

0m 

1.7 

1.704 

0s 

1.3 

1.3001 

2.4.  Numerical  method 


The  governing  equations  are  discretized  according  to  the  finite 
volume  method.  A  Newton-Raphson  solver  which  implements 
the  generalized  minimal  residual  method  (GMRES)  [32]  with  a 
localized  Incomplete  LU  preconditioner  [33]  is  used  to  solve  the 
discretized  governing  equations.  The  simulations  are  parallelized 
using  the  Message  Passing  Interface  (MPI)  Library  [34]  and  are  run 
on  16  processors.  The  simulations  are  considered  to  be  converged 
when  the  nonlinear  rms  is  converged  to  machine  zero.  The  consis¬ 
tency  of  the  simulation  results  are  verified  by  comparing  the  total 
oxygen  consumption,  total  proton  consumption,  total  electron  con¬ 
sumption  and  total  water  production  values.  More  details  on  the 
numerical  approach  for  a  similar  problem  are  given  in  [35]. 

The  effective  transport  parameters  for  the  simulation  are  com¬ 
puted  in  the  following  manner.  The  total  flux  for  any  given  quantity 
g  (heat  or  species)  can  be  represented  by  an  effective  transport 
parameter  Mep  the  length  of  the  solution  domain  l,  and  the  spec¬ 
ified  values  at  opposite  boundaries  g\  and  g2.  Using  this  notation, 
the  total  flux  through  the  solution  domain  Ag  is  expressed  as 

Ag  =  -Meffi82-g'].  (21) 

The  specified  boundary  conditions  and  length  of  the  domain  are 
known.  The  total  flux  through  the  domain  Ag  can  be  computed 
from  the  simulation.  Thus,  the  effective  transport  parameter  can 
be  computed  as 


Meff  =  - 


Agl 


(&  -  gi ) 


(22) 


The  computation  of  the  flux  Ag  is  rather  straightforward  when 
electrochemical  reactions  are  not  considered  in  the  model.  In  this 
case,  the  influx  at  one  boundary  must  equal  to  the  outflux  at  the 
opposite  boundary.  However,  when  electrochemical  reactions  are 
considered  in  the  model,  there  is  net  production  or  consumption  in 
the  domain,  and  the  total  influx  and  outflux  values  are  not  equal. 
In  the  current  work,  it  was  observed  that  the  differences  between 
the  influx  and  outflux  values  for  each  species  are  quite  small  and 
thus  it  is  reasonable  to  use  the  outflux  values  when  computing  the 
effective  transport  parameters. 


3.  Results 


Simulations  are  run  for  the  aforementioned  16  data  sets  for 
a  wide  variety  of  random  microstructures  which  have  pore  vol¬ 
ume  fractions  between  25  and  50  percent.  For  each  reconstructed 
microstructure,  simulation  results  are  computed  for  three  orienta¬ 
tions.  Each  orientation  has  specified  dirichlet  boundary  conditions 
at  different  opposing  boundaries:  one  for  the  x  boundaries  (on  the 
y-z  planes),  one  for  they  boundaries  (on  the  x-z  planes)  and  one  for 
the  Z  boundaries  (on  the  x-y  planes).  The  consumption  values  and 
effective  transport  parameters  computed  from  these  three  simula¬ 
tions  are  averaged.  It  is  the  average  values  that  are  represented  in 
subsequent  plots. 

Given  that  the  operating  conditions  for  each  case  are  different, 
these  correspond  to  different  current  densities  that  would  be  found 
in  a  fuel  cell.  However,  only  a  200  nm  x  200  nm  x  200  nm  section  of 
the  catalyst  layer  is  considered,  and  thus  precise  values  for  fuel  cell 
current  densities  cannot  be  obtained.  In  particular,  at  low  relative 
humidities,  it  is  likely  that  the  electrochemical  reaction  will  only 
take  place  in  a  very  thin  layer  of  the  catalyst  layer  that  is  close  to  the 
membrane.  Once  all  of  the  incoming  protons  have  been  consumed 
through  electrochemical  reactions,  the  remainder  of  the  catalyst 
layer  will  be  inactive. 

However,  one  can  obtain  approximate  current  density  ranges  for 
the  given  cases.  In  this  case,  one  can  convert  the  total  oxygen  con- 
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Table  5 

Relative  humidity  values  and  approximate  current  density  ranges  for  considered 
cases. 


T(  K) 

Ch2o  (mol  nrr3) 

Co2  (mol  nrr3) 

a 

Current  density  (A cm-2) 

343 

10 

14 

0.91 

0.75-1.5 

343 

10 

7 

0.91 

0.38-0.75 

343 

8 

14 

0.73 

0.33-0.66 

343 

8 

7 

0.73 

0.17-0.33 

353 

10 

14 

0.62 

0.33-0.61 

353 

10 

7 

0.62 

0.17-0.30 

353 

8 

14 

0.49 

0.06-0.17 

353 

8 

7 

0.49 

0.03-0.08 

sumption  values  (mol  s-1 )  into  the  total  reaction  current  in  A  cm-3. 
In  order  to  estimate  the  current  density  for  a  given  case,  this  reac¬ 
tion  current  can  be  multiplied  by  the  effective  reaction  thickness  of 
the  catalyst  layer.  The  simulation  results  indicate  that  the  limiting 
factor  in  catalyst  layer  utilization  is  the  low  membrane  conductiv¬ 
ity,  and  this  can  then  be  used  to  determine  an  effective  reaction 
thickness.  Thus,  for  the  case  with  the  highest  relative  humidity 
value  (T  =  343  K,  cH2o  =  1 0  mol  cm-3 ),  it  is  assumed  that  this  length 
is  between  half  the  catalyst  layer  thickness  and  the  entire  catalyst 
layer  thickness  (taken  to  be  10  [xm  here).  Thus,  an  approximate 
current  density  range  can  be  obtained  for  this  case.  For  cases  with 
lower  relative  humidities,  a  factor  is  introduced  based  on  the  ratio 
of  the  membrane  conductivity  to  the  conductivity  at  the  maximum 
relative  humidity. 

Based  on  these  calculations,  Table  5  lists  the  approximate  range 
of  relative  humidities,  conversion  factors,  and  approximate  range  of 
current  densities  for  each  case.  Because  the  membrane  conductivity 
is  computed  based  on  experimental  data  taken  for  a  specific  recast 
membrane,  the  current  densities  at  low  relative  humidity  values  are 
very  small.  These  values  will  change  when  different  membranes  are 
used. 

3.1.  Effective  transport  parameters 

Previous  modeling  work  [13]  indicated  that  when  larger  car¬ 
bon  spheres  are  used,  larger  pores  are  created  between  the  carbon 
spheres  which  results  in  higher  Knudsen  diffusivities  and  hence 
higher  effective  oxygen  diffusivities.  The  same  phenomenon  is 
observed  here,  as  shown  in  Fig.  2.  There  was  very  little  change 
observed  in  the  effective  oxygen  diffusivity  when  the  tempera¬ 
ture  boundary  conditions  were  changed  within  a  normal  operating 
range  (343-353  K). 


Fig.  2.  Effective  diffusivities  of  oxygen  for  cases  with  different  sphere  radii.  The  four 
numbers  in  the  series  label  correspond  to  T (I<),  R  (mol  nrr3),  Psat  (mol  irr3),  sphere 
radius  (nm). 


Fig.  3.  Effective  diffusivities  of  water  vapor  for  cases  with  different  sphere  radii 
and  water  vapor  concentrations.  Labels:  T(K),Ch2o  (mol  rrr3),  co2  (mol  nrr3),  sphere 
radius  (nm). 

Fig.  3  shows  the  computed  effective  water  vapor  diffusivity  for 
four  different  sets  of  data.  Similar  to  the  effective  oxygen  diffu¬ 
sivity  data,  when  larger  carbon  spheres  are  used,  larger  effective 
diffusivities  are  obtained  due  to  larger  pores  and  higher  Knudsen 
numbers.  Another  interesting  feature  of  the  data  is  that  at  low  pore 
volume  fractions,  data  sets  with  higher  water  vapor  concentrations 
have  higher  effective  diffusivities.  This  is  due  to  the  phenomenon 
of  electro-osmotic  drag.  Since  the  potential  difference  across  the 
catalyst  layer  sample  is  prescribed,  at  higher  water  vapor  concen¬ 
trations,  there  is  increased  proton  transport  through  the  membrane 
due  to  higher  membrane  conductivities.  Water  vapor  molecules 
attach  to  the  protons  and  are  carried  in  this  manner.  This  effect 
is  most  pronounced  at  low  pore  volume  fractions,  because  the 
ionomer  volume  fraction  is  higher  and  higher  effective  proton  con¬ 
ductivity  values  are  obtained. 

Effective  proton  conductivity  data  for  a  wide  range  of  cases  is 
shown  in  Fig.  4.  As  expected,  the  most  important  factor  in  deter¬ 
mining  the  effective  proton  conductivity  is  the  relative  humidity, 
which  depends  on  the  temperature  and  the  water  vapor  concen¬ 
tration  as  shown  in  Eq.  (16).  When  larger  spheres  are  used,  for  a 
given  pore  volume  fraction,  the  carbon  volume  fraction  is  higher 
and  the  ionomer  volume  fraction  is  lower,  hence  the  data  for  the 
larger  spheres  corresponds  to  lower  ionomer  volume  fraction  val¬ 
ues.  The  effective  electron  conductivity  data  is  shown  in  Fig.  5.  The 


Fig.  4.  Effective  proton  conductivities  for  cases  with  different  sphere  radii,  water 
vapor  concentrations,  and  temperatures.  Labels :  T ( K),  ch2 o  ( mol  irr 3 ),  co2  ( mol  m_  3 ), 
sphere  radius  (nm),  relative  humidity. 
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Fig.  5.  Effective  electron  conductivities  for  cases  with  different  sphere  radii.  Labels: 
T(K),  ch2o  (molm-3),  Co2  (molm-3),  sphere  radius  (nm). 

effective  electron  conductivities  are  much  larger  than  the  effec¬ 
tive  proton  conductivities,  because  the  bulk  electron  conductivity 
is  much  larger.  There  is  also  much  more  variance  in  the  effective 
electron  conductivity  data.  This  may  be  due  to  the  low  overlap  tol¬ 
erance  and  the  presence  of  isolated,  disconnected  carbon  spheres 
that  are  present  in  some  microstructures. 

Figs.  6  and  7  show  the  effective  thermal  conductivity  as  a 
function  of  pore  volume  fraction  and  ionomer  volume  fraction, 
respectively.  As  the  porosity  of  the  catalyst  layer  is  increased,  the 
thermal  conductivity  decreases,  because  the  thermal  conductivity 
of  air  is  the  lowest  of  the  three  phases.  The  amount  of  ionomer 
in  the  catalyst  layer  increases  with  the  number  of  carbon  spheres, 
and  thus  the  thermal  conductivity  increases  with  an  increase  in  the 
ionomer  volume  fraction. 

3.2.  Consumption  plots 

Figs.  8  and  9  show  the  total  oxygen  consumption  in  the  cata¬ 
lyst  layer  sample  for  cases  run  at  343  and  353  K.  It  is  clear  from 
the  data  that  at  higher  oxygen  concentrations  and  higher  tempera¬ 
tures,  higher  consumption  values  are  obtained.  This  is  consistent 
with  the  Tafel  equation  and  the  expressions  for  exchange  cur¬ 
rent  density  and  the  charge  transfer  coefficient  as  given  in  Eqs. 
(6)  and  (7).  Another  dominant  trend  is  that  as  the  pore  volume 
fraction  increases,  the  total  consumption  decreases.  This  is  due 


Fig.  6.  Effective  thermal  conductivities  as  a  function  of  pore  volume  fraction.  Labels: 
T(K),  Ch2o  (molm-3),  Co2  (molm-3),  sphere  radius  (nm). 


Fig.  7.  Effective  thermal  conductivities  as  a  function  of  ionomer  volume  fraction. 
Labels:  T (K),  ch2o  (molm-3),  co2  (molm-3),  sphere  radius  (nm). 


Fig.  8.  Total  oxygen  consumption  as  a  function  of  pore  volume  fraction  at  343  K. 
Labels:  T (K),  cH2o  (molm-3),  c02  (molm-3),  sphere  radius  (nm). 

to  the  fact  that  there  are  more  carbon  spheres  in  the  microstruc¬ 
tures  at  low  pore  volume  fractions  and  thus  more  area  available  for 
electrochemical  reactions.  When  larger  carbon  spheres  are  used, 
there  is  less  surface  area  for  reactions  at  a  given  carbon  volume 
fraction  than  when  smaller  carbon  spheres  are  used.  Thus,  lower 


Fig.  9.  Total  oxygen  consumption  as  a  function  of  pore  volume  fraction  at  353  K. 
Labels:  T(K),  Ch2o  (molm-3),  Co2  (molm-3),  sphere  radius  (nm). 
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Fig.  10.  Comparison  of  simulation  data  with  experimental  data  [36]  for  effective 
oxygen  diffusivity. 

consumption  values  are  obtained  when  larger  carbon  spheres  are 
used. 

4.  Discussion 

Comparing  the  results  obtained  from  the  computational  model 
with  published  experimental  data  yields  further  insights.  Recently, 
the  effective  oxygen  diffusivity  through  a  PEM  fuel  cell  catalyst 
layer  was  measured  by  researchers  from  General  Motors  Corpo¬ 
ration  [36]  at  a  temperature  of  353  K  and  a  pressure  of  1  atm. 
Their  values  are  much  lower  than  those  found  in  the  simulations, 
as  shown  in  Fig.  10. 

The  trends  are  similar,  but  it  appears  that  the  effective  diffusivity 
curves  found  from  the  simulations  are  shifted  to  the  left  relative  to 
the  experimental  data.  A  number  of  reasons  could  account  for  this 
shift.  Liquid  water  is  not  accounted  for  in  the  model  and  this  may 
present  an  impediment  to  oxygen  diffusion.  Also,  ionomer  swelling 
is  not  taken  into  account  in  the  model.  The  swelling  of  the  ionomer 
film  could  constrict  the  pores,  resulting  in  lower  pore  diameters  and 
thus  lower  effective  diffusivities.  However,  it  is  unlikely  that  these 
aspects  can  alone  explain  the  difference  in  the  results,  because  the 
experimental  results  show  only  a  small  change  in  the  measured 
effective  diffusivities  when  the  relative  humidity  is  changed  from 
0  to  100  percent. 

An  alternative  explanation  for  the  discrepancy  is  that  the 
assumption  of  a  uniformly  thin  ionomer  film  is  not  valid.  A  uni¬ 


formly  thin  layer  around  the  carbon  spheres  does  not  present  a 
significant  resistance  to  oxygen  diffusion  through  the  catalyst  layer, 
since  the  oxygen  already  has  to  travel  around  the  carbon  spheres. 
Some  self-organization  studies  have  predicted  that  the  surface  of 
the  carbon  particles  is  not  completely  covered  by  the  ionomer  [37]. 
If  the  ionomer  is  distributed  nonuniformly,  this  might  present  more 
of  a  barrier  for  oxygen  diffusion  as  ionomer  filaments  and  agglom¬ 
erates  may  protrude  into  pores  between  the  carbon  spheres.  An 
example  of  this  is  shown  in  Fig.  1 1 .  The  nonuniform  ionomer  distri¬ 
bution  would  significantly  change  the  pore  morphology,  increase 
the  tortuosity  of  diffusing  gases  and  decrease  their  effective  dif¬ 
fusivities.  Verification  of  this  conjecture  has  to  await  progress  in 
detailed  characterization  of  the  catalyst  layer  morphology. 

An  additional  assumption  of  the  model  is  that  the  ionomer  dis¬ 
tribution  and  porosity  throughout  the  length  of  the  catalyst  layer 
are  approximately  uniform.  This  assumption  allows  computation 
of  the  effective  transport  properties  from  a  computational  domain 
with  dimensions  that  are  2  orders  of  magnitude  smaller  than  the 
thickness  of  most  catalyst  layers.  If,  however,  there  was  more 
ionomer  distributed  towards  one  end  of  the  catalyst  layer  than 
the  other  (meaning  that  the  porosity  through  the  catalyst  layer 
also  varies),  this  would  represent  a  significant  barrier  for  oxygen 
transport.  The  measured  effective  oxygen  diffusivity  throughout 
the  entire  catalyst  layer  would  be  much  lower  for  a  catalyst  layer 
with  inhomogeneous  porosities  and  ionomer  distributions. 

Finally,  the  Knudsen  diffusivity  is  derived  for  one-dimensional 
transport  in  a  straight  channel.  Every  molecule  that  collides  with 
the  wall  changes  its  direction  but  continues  to  move  towards 
regions  of  lower  concentrations.  In  the  case  of  a  PEM  fuel 
cell  catalyst  layer,  the  microstructure  is  composed  of  randomly 
placed  carbon  spheres.  Molecular  collisions  may  redirect  molecules 
towards  regions  of  higher  concentrations  and  this  effect  would 
significantly  reduce  the  effective  diffusivity.  Furthermore,  oxy¬ 
gen  molecules  in  the  catalyst  layer  “collide”  with  an  ionomer 
layer  which  is  permeable  to  oxygen.  This  is  different  from  classi¬ 
cal  Knudsen  diffusion  where  molecules  collide  with  impermeable 
structures. 

An  interesting  insight  from  these  results  comes  from  the  fact 
that  agglomerates  are  not  explicitly  formed  in  the  reconstruction 
process.  Each  carbon  sphere  is  placed  individually  in  the  domain, 
resulting  in  a  unimodal  pore  size  distribution.  An  alternative 
approach  would  be  to  form  agglomerates  with  a  specified  num¬ 
ber  of  carbon  spheres,  and  form  the  catalyst  layer  microstructure 
by  placing  agglomerates  in  the  domain,  covering  the  agglomerates 
with  ionomer.  This  would  create  a  bimodal  pore  size  distribution  in 
the  domain  due  to  micropores  inside  the  agglomerates  and  macro¬ 
pores  between  the  agglomerates.  However,  if  such  a  scheme  were 


(a)  Uniform  ionomer  distribution  (b)  Nonuniform  ionomer  distribution 


Fig.  11.  Schematic  showing  diffusion  paths  for  oxygen  in  microstructures  with  different  ionomer  distributions.  It  is  clear  that  the  oxygen  in  Fig.  lib  has  a  more  tortuous 
path  to  travel. 
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Fig.  12.  Comparison  of  simulation  data  with  experimental  data  for  the  effective 
proton  conductivity  of  the  catalyst  layer.  Measurements  are  assumed  to  take  place 
at  353  K  unless  otherwise  specified.  The  data  by  Iden  et  al.  was  taken  for  two  pseudo 
catalyst  layers:  one  with  Ketjen  Black  (KB)  as  the  carbon  support  and  one  with 
graphitized  Ketjen  Black  (GKB)  as  the  carbon  support. 


implemented  in  the  current  model,  the  effective  diffusivities  of  oxy¬ 
gen  would  be  higher  due  to  the  presence  of  macropores,  which 
would  reduce  the  impact  of  Knudsen  diffusion.  Thus,  for  a  model 
which  considered  agglomerates,  there  would  be  an  even  larger  dif¬ 
ference  between  experimental  and  computational  results. 

A  number  of  experimental  works  have  sought  to  measure  the 
effective  proton  conductivity  or  resistance  through  the  catalyst 
layer  of  a  PEM  fuel  cell  [38-42].  In  general,  the  experimental 
data  has  shown  a  linear  relationship  between  the  effective  pro¬ 
ton  conductivity  and  the  ionomer  volume  fraction  [38,41],  which 
is  qualitatively  similar  to  the  computational  results,  as  shown  in 
Fig.  12.  However,  the  normalized  effective  conductivity  values  com¬ 
puted  by  the  simulations  for  a  given  ionomer  volume  fraction 
are  much  lower  than  values  obtained  from  experiments,  with  the 
exception  of  the  data  taken  from  the  catalyst  with  graphitized  Ket¬ 
jen  Black  (GKB). 

The  discrepancy  between  experimental  and  computational 
results  again  calls  into  question  the  aforementioned  assumption  of 
uniformly  thin  ionomer  layers.  This  assumption  requires  the  pro¬ 
tons  to  travel  close  to  the  surface  of  the  carbon  spheres  and  the 
tortuosity  may  be  decreased  if  there  are  more  direct  routes  that 
protons  can  take.  These  may  be  created  by  having  thicker  ionomer 
regions  at  the  intersections  between  two  carbon  spheres.  Another 
consideration  is  that  the  current  model  uses  an  approximately  con¬ 
stant  ionomer/carbon  (I/C)  ratio  for  a  given  carbon  sphere  radius 
and  ionomer  thickness  regardless  of  the  porosity  (1.1  for  carbon 
sphere  radii  of  1 2  nm  and  0.7  for  carbon  sphere  radii  of  20  nm  with 
ionomer  thicknesses  of  4  nm).  The  model  also  assumes  that  the  car¬ 
bon  spheres  all  have  the  same  radii.  It  is  more  likely  that  there  is 
a  distribution  with  regard  to  the  size  of  the  carbon  particles,  and 
in  a  previous  work,  it  was  observed  that  the  ionomer  conductivity 
increased  when  a  normal  distribution  of  carbon  particles  was  used 
[13]. 

However,  the  fact  that  the  experimental  data  yields  tortuosity 
factors  that  are  less  than  or  close  to  one  indicates  that  there  must  be 
a  different  mechanism  of  proton  transport  in  thin  ionomer  films  on 
carbon  than  what  is  measured  with  bulk  conductivity  values.  Tortu¬ 
osity  values  of  less  than  one  are  unphysical.  It  is  also  highly  unlikely 
that  the  actual  tortuosity  is  close  to  one  because  this  implies  that 
there  are  straight  proton-conducting  ionomer  “tubes”  which  span 
the  entire  thickness  of  the  catalyst  layer  from  one  side  to  the  other. 


There  are  several  possibilities  as  to  what  this  different  mech¬ 
anism  might  be,  with  interactions  between  the  carbon  support 
material  and  the  ionomer  being  one  possibility.  The  proton  con¬ 
ductivity  of  the  ionomer  has  been  shown  to  depend  on  the  carbon 
material  [42]  and  on  the  I/C  ratio  [41].  Indeed,  results  from  the 
computational  model  are  very  close  to  experimental  data  obtained 
using  GKB  as  a  carbon  support  [42].  An  additional  possibility  is  that 
some  mechanism  of  surface  diffusion  takes  place  as  there  is  a  much 
higher  surface  to  volume  ration  for  thin  ionomer  films  than  for  bulk 
films.  Finally,  some  have  suggested  that  proton  transport  in  the 
catalyst  layer  occurs  through  water  [37]  and  this  additional  trans¬ 
port  mechanism  may  serve  to  increase  the  total  effective  proton 
conductivity  through  the  catalyst  layer. 

The  insights  that  have  been  obtained  from  comparing  the  results 
of  this  model  with  experimental  results  can  be  used  to  guide  the 
future  development  of  catalyst  layer  models  and  raise  a  number 
of  questions.  Research  must  be  done  to  obtain  an  expression  for 
the  Knudsen  diffusivity  which  is  more  appropriate  for  catalyst 
layer  microstructures.  In  addition,  different  reconstruction  meth¬ 
ods  must  be  investigated  to  determine  the  effect  of  the  ionomer 
distribution  on  the  effective  oxygen  diffusivity.  To  improve  the 
model  of  proton  transport,  experimental  proton  conductivity  val¬ 
ues  should  be  obtained  for  thin  ionomer  films  on  carbon  supports 
instead  of  bulk  ionomer  membranes. 

Future  development  of  the  catalyst  layer  model  will  incorporate 
the  discrete  placement  of  platinum  particles  and  a  correlation  with 
the  platinum  loading,  ionomer  swelling,  and  the  presence  of  liquid 
water.  In  addition,  the  model  will  account  for  multispecies  diffusion 
and  for  variations  in  the  thermal  conductivity  of  the  ionomer  due 
to  changes  in  the  temperature  and  relative  humidity. 

5.  Conclusions 

A  fully  coupled  model  has  been  developed  for  the  simulation  of 
oxygen,  protons,  electrons,  water  vapor  and  temperature  in  a  small 
section  of  the  catalyst  layer  of  a  PEM  fuel  cell.  Effective  transport 
parameters  and  consumption  values  are  computed  for  a  large  range 
of  different  microstructures  for  a  wide  variety  of  operating  condi¬ 
tions.  The  differences  between  the  simulation  results  and  published 
experimental  results  are  discussed  and  the  model  assumptions 
are  critically  re-examined  to  provide  insight  and  guidance  for  the 
development  of  accurate  catalyst  layer  simulations. 
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Glossary 


li„t:  term  is  equal  to  one  at  carbon/ionomer  interface  and  zero  elsewhere 
a:  relative  humidity 

Cj :  coefficient  for  membrane  conductivity  (S  cm-1 ) 
c2:  coefficient  for  membrane  conductivity  (I<-1 ) 


C3:  coefficient  for  membrane  conductivity  (K-2) 

C4:  coefficient  for  membrane  conductivity  (K-3) 

C5:  coefficient  for  membrane  conductivity  (K-4) 

Cg.'  coefficient  for  membrane  conductivity  (S  cm-1 ) 

Ch2o •'  water  vapor  concentration  (mol cm-3) 

Co2:  oxygen  concentration  (mol cm-3) 

Co2,ref  reference  oxygen  concentration  (mol cm-3) 
dp:  pore  diameter  (cm) 

Dh2o-  water  vapor  diffusivity  (cm2  s-1 ) 

%2o,b;  bulk  water  vapor  diffusivity  in  pore  region  (cm2  s_1 ) 
DH2o,i •'  ionomer  water  vapor  diffusivity  (cm2  s-1 ) 

Dh2o,kti :  Knudsen  water  vapor  diffusivity  in  pore  region  (cm2  s-1 ) 
Do2  :  oxygen  diffusivity  (cm2  s_1 ) 

D02y.  bulk  oxygen  diffusivity  in  pore  region  (cm2  s_1 ) 

Do2)j.'  ionomer  oxygen  diffusivity  (cm2  s-1 ) 

Do2,Kn:  Knudsen  oxygen  diffusivity  in  pore  region  (cm2  s-1 ) 

F:  Faraday’s  constant  (C  mol-1 ) 

g:  solution  variable  for  quantity  in  the  domain 

gj:  solution  variable  at  boundary  1 

g2:  solution  variable  at  boundary  2 

i0:  exchange  current  density  (A cm-2) 

l:  length  of  the  solution  domain  (cm) 

Meff:  effective  transport  parameter 
nd:  electro-osmotic  drag  coefficient 
p;  pressure  (atm) 

p2:  coefficient  for  saturation  pressure  (Pa) 
p2:  coefficient  for  saturation  pressure  (Pa  K_1 ) 
p3:  coefficient  for  saturation  pressure  (Pa  K-1 ) 
p4:  coefficient  for  saturation  pressure  (Pa  K-1 ) 

Psat:  saturation  pressure  (Pa) 

R:  universal  gas  constant  (J  mol-1  K_1 ) 

AS/,:  change  in  entropy  (J  mol-1  K-1 ) 

T:  temperature  (K) 

Greek  symbols 

ac:  charge  transfer  coefficient 
p:  overpotential  (V) 

X:  thermal  conductivity  (W  cm-1  K-1 ) 

Ag:  total  flux  through  the  domain 
77:  Peltier  coefficient  (J  C_1 ) 
am:  membrane  conductivity  (S  cm-1 ) 
os:  solid  conductivity  (S  cm-1 ) 

(f)m:  membrane  potential  (V) 

(j)s:  solid  potential  (V) 


